function [allMSPs] = calculate_OPALM_MASTER(data,pressures,TemperatureIndicies,CompIndicies,ModelVariablesIndicies, garnet_coefficients, spinel_coefficients, plagioclase_coefficients)
%The process is to:
%(1) first calculate the pressure
%(2) then using that P calculate T and X
%(3) calculate the X difference to see if the liquid is is EQ with a
%lherzolite



%%%%%%%%%%%%%%%%%%%%%%%%
% [A,TemperatureIndicies] = ismember(VariablesT, ElementsZZZ);
% [A,CompIndicies] = ismember(VariablesX, ElementsZZZ);
% 
% 
% FittedTargetVariables = {'P' 'T' 'Qtz' 'Plag' 'Oliv' 'Cpx'};
% ModelVariables = {'P' 'T' 'Oliv' 'Cpx' 'Plag' 'Qtz' 'P'};
% [A,ModelVariablesIndicies] = ismember(FittedTargetVariables, ModelVariables);




data(isnan(data)) = 0;

for i = 1:size(data,1)
    
    %variableValues = [onesmatrix data(:,PressureIndicies)]; %add in a ones column account for the intercept
    
    plagioclasePressures = pressures;
    
    
    % garnetPressures = deepP:-1:spinel2garnet(i);
    % spinelPressures = spinel2garnet(i):-1:plag2spinel(i);
    % plagioclasePressures = plag2spinel(i):-1:shallowP;
    
    
    % %Step 2 : calculate T and X
    onesmatrix = ones(size(plagioclasePressures,2),1);
    % variableValues = [onesmatrix (plagioclasePressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
    % plagioclaseTX= [plagioclasePressures' variableValues*plagioclase_coefficients'];
    
    %Step 2 : calculate T and X
    variableValues = [onesmatrix (plagioclasePressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
    T = [plagioclasePressures' variableValues*plagioclase_coefficients'];
    
    variableValues = [onesmatrix (plagioclasePressures)' repmat(data(i,CompIndicies),size(onesmatrix,1),1)];
    X = [plagioclasePressures' variableValues*plagioclase_coefficients'];
    
    combinedTX=[T(:,1:2), X(:,3:end)];
    
    singleDataPointResults = [combinedTX];
    
    
    resultsgood = singleDataPointResults(:,ModelVariablesIndicies);
    if isempty(resultsgood == 1)
        resultsgood =  NaN.*zeros(41,6);
    end
    
    
    allMSPs(:,:,i) =resultsgood;
    
    
    
    
end
end

